# This file produces tables 5 and 6
set.seed(123)
setwd("/Users/ericguntermann/Documents/Papers/Representation of Party Preferences/Replication")

## Criterion 2 (gallagher)
library(foreign)
criteria <- read.dta("criteria.dta")
d1 <- with(criteria, list(mostliked=mostliked, gallagher=log(gallagher), gdppercap=gdppercap, freedomhouse=freedomhouse))
d1$N <- length(unique(criteria$countryyear))
m1 <- function(){
  for(i in 1:N){
    mostliked[i] ~ dbern(p[i]) 
    probit(p[i]) <- mu[i]
    mu[i] <- b + b.gallagher*gallagher[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i]
  }
  b ~ dnorm(0,0.01)
  b.gallagher ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
}

i1 <- function(){
  list("b"=rnorm(1), "b.gallagher"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1))
}
p1 <- c("b", "b.gallagher", "b.gdppercap", "b.freedomhouse")


## Criterion 2 (mdm)
library(foreign)
criteria <- read.dta("criteria.dta")
d2 <- with(criteria, list(mostliked=mostliked, mdm=log(mdm), gdppercap=gdppercap, freedomhouse=freedomhouse))
d2$N <- length(unique(criteria$countryyear))

m2 <- function(){
  for(i in 1:N){
    mostliked[i] ~ dbern(p[i]) 
    probit(p[i]) <- mu[i]
    mu[i] <- b + b.mdm*mdm[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i]
  }
  b ~ dnorm(0,0.01)
  b.mdm ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
  
}

p2 <- c("b", "b.mdm", "b.gdppercap", "b.freedomhouse")


i2 <- function(){
  list("b"=rnorm(1), "b.mdm"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1))
}



## Criterion 2 (gallagher) with number of parties
library(foreign)
dat <- read.dta("criteria.dta")
d3 <- with(dat, list(mostliked=mostliked, gallagher=log(gallagher), gdppercap=gdppercap, freedomhouse=freedomhouse, ngov=ngov))
d3$N <- length(unique(dat$countryyear))
m3 <- function(){
  for(i in 1:N){
    mostliked[i] ~ dbern(p[i]) 
    probit(p[i]) <- mu[i]
    mu[i] <- b + b.gallagher*gallagher[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i] + b.ngov*ngov[i]
  }
  b ~ dnorm(0,0.01)
  b.gallagher ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  b.ngov ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
}

i3 <- function(){
  list("b"=rnorm(1), "b.gallagher"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1), "b.ngov"=rnorm(1))
}
p3 <- c("b", "b.gallagher", "b.gdppercap", "b.freedomhouse", "b.ngov")



## Criterion 2 (mdm) with number of parties
library(foreign)
criteria <- read.dta("criteria.dta")
d4 <- with(criteria, list(mostliked=mostliked, mdm=log(mdm), gdppercap=gdppercap, freedomhouse=freedomhouse, ngov=ngov))
d4$N <- length(unique(criteria$countryyear))
m4 <- function(){
  for(i in 1:N){
    mostliked[i] ~ dbern(p[i]) 
    probit(p[i]) <- mu[i]
    mu[i] <- b + b.mdm*mdm[i] + b.gdppercap*gdppercap[i] + b.freedomhouse*freedomhouse[i] + b.ngov*ngov[i]
  }
  b ~ dnorm(0,0.01)
  b.mdm ~ dnorm(0,0.01)
  b.gdppercap ~ dnorm(0,0.01)
  b.freedomhouse ~ dnorm(0,0.01)
  b.ngov ~ dnorm(0,0.01)
  tau ~ dgamma(1,1)
  
}

p4 <- c("b", "b.mdm", "b.gdppercap", "b.freedomhouse", "b.ngov")


i4 <- function(){
  list("b"=rnorm(1), "b.mdm"=rnorm(1), "b.gdppercap"=rnorm(1), "b.freedomhouse"=rnorm(1), "b.ngov"=rnorm(1))
}





d <- list(d1,d2,d3,d4)
p <- list(p1,p2,p3,p4)
m <- list(m1,m2,m3,m4)
ini <- list(i1,i2,i3,i4)


library(foreach)
library(doParallel)
library(R2jags)
library(doRNG)
registerDoParallel(4) 
cl <- makeCluster(4)
registerDoParallel(cl)
results <- foreach(i=1:4, .packages = c('R2jags')) %dopar% jags(data=d[[i]], model=m[[i]], parameters.to.save=p[[i]], inits=ini[[i]], n.iter=100000, n.burnin=50000, n.thin=10, n.chains=2)


results1 <- as.mcmc(results[[1]])
results2 <- as.mcmc(results[[2]])
results3 <- as.mcmc(results[[3]])
results4 <- as.mcmc(results[[4]])

# Convergence diagnostics

library(ggmcmc)
p1 <- ggs_geweke(ggs(results1)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 1")
p2 <- ggs_geweke(ggs(results2)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 2")
p3 <- ggs_geweke(ggs(results3)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 3")
p4 <- ggs_geweke(ggs(results4)) + theme(legend.position = "none") + ylab("") + ggtitle("Model 4")
library(gridExtra)

pdf("criterion2_conv.pdf")
grid.arrange(p1,p2,p3,p4)
dev.off()

## GET OUTPUT


#Model 1 (Gallagher)
coefs1 <- as.data.frame(as.matrix(results1))
b <- coefs1[,"b"]
b.gallagher <- coefs1[,"b.gallagher"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]


co1 <- c(mean(b), mean(b.gallagher), NA, mean(b.freedomhouse), mean(b.gdppercap))
sd1 <- c(sd(b),sd(b.gallagher), NA, sd(b.freedomhouse), sd(b.gdppercap))
p1 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), ifelse(mean(b.gallagher)>0, mean(b.gallagher>0),mean(b.gallagher<0)), NA, ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))


#Model 2 (MDM)
coefs1 <- as.data.frame(as.matrix(results2))
b <- coefs1[,"b"]
b.mdm <- coefs1[,"b.mdm"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]

co2 <- c(mean(b), NA, mean(b.mdm), mean(b.freedomhouse), mean(b.gdppercap))
sd2 <- c(sd(b), NA, sd(b.mdm), sd(b.freedomhouse), sd(b.gdppercap))
p2 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), NA, ifelse(mean(b.mdm)>0, mean(b.mdm>0),mean(b.mdm<0)), ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))


# Table 5
library(xtable)
model.effects <- data.frame(co1,sd1,p1,co2,sd2,p2)                                                                                                                
rownames(model.effects) <- c("Intercept", "Gallagher*", "MDM*", "Freedom House","GDP per capita")
colnames(model.effects) <- c("Mean", "SD", "P", "Mean", "SD", "P")
xtab <- xtable(model.effects, caption="Criterion 2", digits=c(2))
sink("criterion2.tex")
print(xtab,caption.placement = getOption("xtable.caption.placement", "top"), table.placement = getOption("xtable.table.placement", "H"), hline.after=0, only.contents=TRUE)
sink(NULL)

#Model 3 (Gallagher with ngov)
coefs1 <- as.data.frame(as.matrix(results3))
b <- coefs1[,"b"]
b.gallagher <- coefs1[,"b.gallagher"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]
b.ngov <- coefs1[,"b.ngov"]


co2 <- c(mean(b), mean(b.gallagher), NA, mean(b.ngov), mean(b.freedomhouse), mean(b.gdppercap))
sd2 <- c(sd(b), sd(b.gallagher), NA, sd(b.ngov),  sd(b.freedomhouse), sd(b.gdppercap))
p2 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), ifelse(mean(b.gallagher)>0, mean(b.gallagher>0),mean(b.gallagher<0)), NA, ifelse(mean(b.ngov)>0, mean(b.ngov>0), mean(b.ngov<0)), ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))


#Model 4 (MDM with ngov)
coefs1 <- as.data.frame(as.matrix(results4))
b <- coefs1[,"b"]
b.mdm <- coefs1[,"b.mdm"]
b.freedomhouse <- coefs1[,"b.freedomhouse"]
b.gdppercap <- coefs1[,"b.gdppercap"]
b.ngov <- coefs1[,"b.ngov"]


co3 <- c(mean(b), NA, mean(b.mdm), mean(b.ngov), mean(b.freedomhouse), mean(b.gdppercap))
sd3 <- c(sd(b), NA, sd(b.mdm), sd(b.ngov), sd(b.freedomhouse), sd(b.gdppercap))
p3 <-  c(ifelse(sd(b)>0, mean(b>0), mean(b<0)), NA, ifelse(mean(b.mdm)>0, mean(b.mdm>0),mean(b.mdm<0)), ifelse(mean(b.ngov)>0, mean(b.ngov>0), mean(b.ngov<0)),ifelse(mean(b.freedomhouse)>0, mean(b.freedomhouse>0), mean(b.freedomhouse<0)), ifelse(mean(b.gdppercap)>0, mean(b.gdppercap>0), mean(b.gdppercap<0)))

# Table 6
library(xtable)
model.effects <- data.frame(co2,sd2,p2,co3,sd3,p3)                                                                                                                
rownames(model.effects) <- c("Intercept", "Gallagher*", "MDM*", "Parties in Govt", "Freedom House","GDP per capita")
colnames(model.effects) <- c("Mean", "SD", "P", "Mean", "SD", "P")
xtab <- xtable(model.effects, caption="Criterion 2", digits=c(2))
sink("criterion2b.tex")
print(xtab,caption.placement = getOption("xtable.caption.placement", "top"), table.placement = getOption("xtable.table.placement", "H"), hline.after=0, only.contents=TRUE)
sink(NULL)



